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ABSTRACT 


An analytic study is made of the flow of a rarefied monatomic 
gas through a two dimensional slot. The parameters of the problem 
are the ratios of downstream to upstream pressures, a, the 
Knudsen number at the high pressure end (based on slot half width) 
Kn Q , and the length to slot half width ratio, Z. First, a moment 
method of solution is used by assuming a discontinuous distribution 
function consisting of four Maxwellians split equally in angular 
space. Numerical solutions are obtained for the resulting equa- 
tions. The characteristics of the transition regime are portrayed 
very well; however the solutions in the free molecule limit are 
systematically lower than the results obtained in that limit by 
more accurate numerical methods. 

Finally, the discrete velocity ordinate method of solution 
is used. The continuous velocity space is represented by 16 
ordinates. Numerical calculations are used to obtain the charac- 
teristics of the transition regime. These characteristics are 
very well represented by this method. Further, the free molecule 
and the slip regime results appear to be highly accurate and serve 
to bolster further confidence in the accuracy of the transition 
regime results. 

For the range of parameters considered, the mass flux through 
the slot, m, is given by. 


m/p 0 a(2RT 0 ) 1/2 


1 r 0.1331 
l 1 "KrT 


Kn. 


(l+a)T-j + 2(T-j - TT 1og s l> 


8/jt 


log s-j log ( 


% + 


A? 


+ 4 


i 



where , 


t _ 0-g) l 

*1 “ l + 4.55 - 2.85 exp{0.2l(Kn o -5)} 

2Kn Q + 8(l+a+T) 

S 1 = 2Kn + 8(l+a-T) 

R is the gas constant, 'a' the slot half width, p Q and T q 
the upstream density and temperature, respectively. The accuracy 
of this formula, compared to the numerical results obtained, is 
5%. It is valid for A = 1,2, 4, 8 and 12, a = 0.1, 0.5 and 0.8 
and Kn Q = °°,5,1 and 0.5. 
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1.0 INTRODUCTION 

From a fluid mechanics viewpoint one of the most attractive 
features of kinetic thoery (with the Boltzmann equation as the 
fundamental equation) is its ability to predict the macroscopic 
properties of a dilute gas in terms of its microscopic properties 
under highly nonequilibrium situations. If 'appropriate' boundary 
conditions are carefully specified, it can be used for both inter- 
nal and external flows. The investigation reported in this thesis 
is a kinetic theory study of the flow of a dilute monatomic gas 
down a finite length (length of slot along direction of flow) slot 
for arbitrary ratios of pressures applied at the ends of the slot. 

The flow first came into focus in the design of pressure 
probes to measure the permeability of lunar surface soil. The 
design features involved the pumping of a gas from a high pressure 
into the soil and carried the assumption that the pressure in the 
flow dropped off to zero far from the probe. If the porous medium 
were visualized as an assembly of cracks and holes, our solution 
could be relevantly applied to the former geometry. It is also 
anticipated that such a flow situation could occur on the Earth's 
soil in places having high pressure gas pockets which are connected 
to the atmosphere through tiny cracks. 

As the mean free path at the high pressure end of the slot is 
varied, it is expected that, if the length of the slot is suffi- 
ciently long, the flow will pass through all stages (or some stages) 
of rarefaction (characterized by the Knudsen number which is the 

ratio of the mean free path of the gas to the slot half width ). 

* 

Half the distance between the slot walls. 
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The above system represents a fundamental flow problem. Rather 
loosely, it can be compared to the problem of the free jet. Whereas 
no solid boundaries are involved there, the slot problem, repre- 
senting an equally nonequilibrium system, has solid boundaries 
playing a decisive role in the flow. The flow is a highly nonlinear 
one, and the number of such highly nonlinear problems that have been 
considered using kinetic theory has been relatively few. In view 
of these considerations, a thorough kinetic theory investigation 
of the flow is considered. 

1.1 Review of Previous Theoretical Work 

In spite of the fundamental nature of the geometry and the 
flow, a solution of sufficient rigor and generality does not seem 
available. 

For free molecule flow through slots of finite length, 
Reynolds and Richey (1967) perform accurate numerical calculations 
for various flow properties. Perhaps the most striking result here 
is that the number flux from the wall varies essentially linearly 
along the length of the slot. A variational solution is also avail- 
able due to Pao and Willis (1962). For flow in the slip regime, 
the only prevailing result is the Poiseuille formula with slip 
boundary conditions. This result, however, is valid only for long 
slots and no corresponding result is available for slots of arbi- 
trary length. 

For slots of infinite length, various calculations in 
the transition regime have been performed. Cercignani and Daneri 
(1963) numerically solved the BGK equation and demonstrated the 
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existence of a 'Knudsen minimum 1 for the linearized flow (with 
small pressure gradients along the slot axis). The same problem 
was also treated by Liu (1968), using the Moment method with a 
choice of the form of the distribution function. It consisted 
of a discontinuous distribution made up of four Maxwellians split 
equally in angular space. Finally, an ad hoe interpolation 
formula, propounded by Hi 1 by and Pahl (1952), is available. It 
involves the interpolation of free molecule and slip flow results 
to transition flow. No comparable results exist for slots of 
finite length. 

1 . 2 Review of Available Theoretical Tools 

The Boltzmann equation is valid for arbitrary deviations 
from equilibrium. However, the complicated mathematical structure 
of the collision term has limited direct use to asymptotic values 
of the Knudsen number. In the last few years a series of different 
methods were used to obtain approximate solutions in the transition 
regime. Some of the methods that could be applied to our problem 
are reviewed briefly below. 

One approach is to use the Moment method. Here the 
Maxwell transport equation replaces the Boltzmann equation as the 
governing equation with the understanding that the exact solution of 
the Boltzmann equation is one that satisfies the Maxwell equation 
for all its moments. However, the solution of an infinite set of 
equations is not feasible. Hence a truncation of the equations is 
performed to get a closed set of equations. This can be done 
through realistic physical assumption--see for example Hamel and 
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Willis (1966). The alternative approach is to recognize that the 
distribution function is formally a function of two sets of inde- 
pendent variables, the physical space position vector and the 
velocity space position vector £j. An explicit form of the 
dependence is assumed, which is compatible with the boundary con- 
ditions. The J*i dependence of the distribution function is still 
an unknown. This is emphasized by the presence of a number of 
unknown functions of in the form of distribution function 
assumed. The number of such unknown functions of FLj required 
should be enough to form a determinate set of equations made up of 
the conservation equations and a reasonable number of higher moment 
equations to represent rarefaction effects. A satisfactory solu- 
tion is hard to define exactly but should certainly represent 
correctly the asymptotic cases of high and low Knudsen number. 

This method is advocated by Lees. A review of the method and appli- 
cations to some flow problems is given in Lees (1965). 

A frequent practice is to replace the collision operator 
in the Boltzmann equation by a relaxation model to give the Bhatnagar, 
Gross and Krook (BGK) equation (1954). Direct numerical solution 
of the BGK equation has been successfully attempted--see for example 
Anderson (1967) and Liepmann, Narasimha and Chahine (1962). For 
linearized internal flows, direct integration of the BGK equation 
has been performed by Cercignani and Daneri (1963). For such flows, 
variational solutions have also been successfully obtained by 
Cercignani and co-workers--see for example Cercignani and Pagani 
(1966). 
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A powerful approach first used in applications in radi- 
ative transfer is the 'discrete ordinate' method. Ithas been 
widely applied by Huang and his co-workers in transition flow, 
both for one and two dimensional, internal and external flow 
problems— see for example Huang and Hartley (1969) and the refer- 
ences cited there. In this method, the continuous velocity space 
is replaced by a number of discrete points to generate a system 
of equations for the velocity distribution function at these dis- 
crete values of the velocity. From the evaluation of these 
distribution functions, any physical quantity can be calculated. 

Also worthy of attention is the 'Restricted Variational 
Method' propounded first by Rosen (1954). Here the task of solving 
the Boltzmann equation is replaced by one of optimizing an integral. 
A form of the distribution function has to be assumed and succeed- 
ing steps closely follow the Rayleigh-Ritz technique. The method 
has since been successfully used by Ortloff (1968) as well. 

A powerful technique for solving the Boltzmann equation 
is the Monte Carlo method. This has been most effectively. employed 
by Bird (1969) and his co-workers. However, the development of 
such a method involves a considerable amount of computer program 
development and comparatively large running times to obtain accu- 
rate results. 

1 . 3 The Present Investigation ■ . . • . 

At the outset a simple discontinuous distribution func- 
tion, made up of four Maxwell ians split equally in angular space, 
is assumed and the Moment method used. The resulting set of semi- 
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linear hyperbolic partial differential equations is integrated 
using the method of characteristics. In the free molecule limit 
(where the Knudsen number based on the slot half width ■* «>) an 
analytical solution is obtained. While the gross features of the 
flow are correctly represented, the free molecule mass flux is 
systematically lower (for the ranges of the various parameters 
considered) than that reported by Reynolds and Richey (1967). 

In fact, for a slot length to half width ratio of 12, the differ- 
ence is almost 35% and it is expected that as this slot length to 
half width ratio increases, the difference would get larger. Also 
the method, as evidenced by the results, seems suitable only for 
slot lengths to half width ratios between 6 and 12. Results in 
the transition regime, for upstream Knudsen numbers less than, or 
equal to, 5 and slot lengths to half width ratios between 1 and 12 
seem reasonable. Also for low values of the Knudsen number, i.e., 
approaching the continuum regime, the results agree remarkably 
well with the Poiseuille flow formula using slip boundary con- 
ditions. 

Due to the deficiency of the Moment method in the free 
molecule limit, a discrete velocity method solution of the BGK 
equation is attempted next. The unknown is a modified distribution 
function. The temperature is set equal to a constant, with negli- 
gible errors expected based on the moment method estimates. The 
velocity space is represented by 16 discrete velocities and the 
distribution function is evaluated at these points by numerical 
calculations. Two distinct methods are used for integrating for 
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the normal mass flux at the wall and the macroscopic quantities 
at any point. The first scheme is a simple trapezoidal rule. The 
second integration scheme allows for the free molecule discontinuity 
of the distribution function and accounts for a linear variation 
along the slot length of the number flux from the walls in the 
free molecule limit. In the free molecule limit, the former method 
leads to overestimates of the mass flux for large slot lengths— as 
much as 35% for a slot length to half width ratio of 12. The latter 
method of integration has errors of less than 0.1% for a slot length 
to half width ratio ranging from 1 to 12. The differences between 
the two modes of integration decrease as the Knudsen number 
decreases, being only 5.0% for length to slot half width ratios 
of 1 to 12, pressure ratios ranging from 0.1 to 0.8, and upstream 
Knudsen number less than 5. 
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2.0 STATEMENT OF THE PROBLEM AND ASSUMPTIONS 

Consider a slot such as is shown pictorially in Fig. la. The 
dimension of the slot la in the direction of the flow is called 
the length of the slot. The dimension of the slot 2a between 
the two walls is called the slot width. The slot has infinite 
depth (the dimension normal to the slot width and length). The 
coordinate system is oriented as shown in Fig. la with walls of 
the slot being represented by yl = ±a, 0 < xl < la. Figure lb 

shows a section (parallel to the flow) of the slot. The slot 
separates two reservoirs containing the same monatomic gas at a 
temperature T . However, on the side xl = 0, the density is 
P 0 » the pressure is p Q ( = p Q RT 0 , with R the gas constant) 
and the Knudsen number is Kn Q (based on the slot half width 'a'). 
On the side xl = al, the pressure is ap Q where a is a non- 
dimensional constant with a value less than 1. 

We make the assumption that the incoming stream on either 
side of the slot is known. These are taken as the distribution 
function in the corresponding reservoirs, i.e., Maxwell ians. This 
is an approximation since the distribution functions of the mole- 
cules coming into the ends of the slot from either reservoir is 
going to be altered by molecules coming out of the slot into the 
reservoirs. There is, however, no easy way of estimating the degree 
of this effect (and hence the errors involved by our assumption) 
unless we solve for the whole flow field. This could be a topic 
for future research. 

Further, there is no accumulation or ablation at the slot 
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walls. Finally, the molecules are assumed to be reemitted from 
the wall diffusely at temperature T . ' 

The problem as posed is to determine the flow field in 
the slot and in particular the mass flux through the slot. Only 
half (above the slot axis yl = 0) of the flow field heed be 
considered due to the symmetry of the flow about the plane yl = 0. 
2.1 Governing Equations and Boundary Conditions 

Let us denote the physical position vector by = 
(xl,yl,zl) and the velocity vector of a molecule by = 

(S xl »5 yl »£ zl ) » where € xl , 5 yl and are the velocity com- 

ponents along the xl , yl and zl axis, respectively. The molecular 
velocity distribution function F, (JR-| ) is defined to be the 

number of molecules per unit of volume of the (FLj ,C.-| ) space. 

The mass density p, velocity U; = (u,v,w), temperature T and 

P.., the flux of * i ' component momentum in the ' j 1 direction, 

^ J 

for the monatomic gas medium, are defined by the following moments 

of F, (&,,£,): 

P(R] ) = m i F-| , 

PM-,) - m /// d 3 ^ £, F 1 

o ? (2J) 

3pRT (jl-j ) = m /// d^ d-, - Ur F ] 

and 

m /// ci 3 ^ 1 ? i Ft , (i = xl ,yl ,zl ) , 

where 

d V= d 5xr<V d5 zi’ ^ 
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m is the mass of a molecule, R is the gas constant and the inte- 
gration is over the whole velocity space. The flow being two 
dimensional, w = 0 and F 1 and hence all its moments are 
functions independent of zl (the coordinate normal to. the xl 
and yl axes). In the steady state with no external force the 
governing equation for can be written as 


Wi 


«F, 

< sr> 


coll 


» 


where V is the spatial gradient operator, and (6F^/6t) C0 -|i is 
the contribution from molecular collisions. In general (6F-j /6t) co -j -| 
is a nonlinear function of F^ . For a dilute gas the Maxwell 
Boltzmann collision integral is the appropriate choice for 
(6F-J / <St) CQ i i (for example, see Kennard, 1938). For mathematical 

simplicity, however, the Boltzmann collision integral is replaced 
by a relaxation type model —the BGK equation (1954). In this model 
equation we write 


6F 1 

coll 0 1 


where v is the collision frequency which is dependent only on 
and 

l“3/2 r tr- ,,\2 


F o = £ (2irRT)“ J/ * exp{-(£ r ur/2RT 0 } 


(2.2a) 


Thus the equation considered is 


3F-| 3F-| 

? xl 9 xT + ^yl 8yT = v[F o' F l ] 


(2.2b) 


In view of the assumptions, the boundary conditions are: 
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At xl = 0, -a < yl * a and £ x i >0, 

F 1 = W~ (2^RT 0 )' 3/2 exp(-^ 2 /2RT 0 ), where ^ 2 = 
At xl = i,*a, -a < yl < a and C xl > °» 

F 1 = IT- ( 2itRT 0 )’ 3/2 exp(-^ 2 /2RT 0 ) ; 

At yl = a, 0 < xl < &*a. 
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(2.2c) 


v = 0, i.e., /// d 3 ^ ? yl F 1 - 0 

and 

F 1 = m ( 2TrRT 0 )" 3/2 e xp(-?i 2 /2RT 0 ) for 5 y] < 0 


At yl = 0, by symmetry about the plane of the slot axis, 
F l^xrSl* 5 zr» x1 «°) - F l^xl^yl^l» x1 ^ 


The collision frequency v is defined as v = pRT/y, 
where the viscosity y is assumed to have the temperature depen- 
dence y/y Q = (T/T 0 ) W , and w is a constants The upstream 
reservoir viscosity y Q is related to the Knudsen number there 
(based on 1 a 1 ) , Kn Q , by 


y o “ 1 p o a K n 0 (8 R T 0 /TT) 1/2 
giving 


V = 


12 RT ° )1/2 < T )1- ( p ) 

2a K "o V p o 


It follows that the local Knudsen number, Kn, 


is given by 


- ( 


f r 1/2 ( 5a ) 


Kn 
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Instead of a Cartesian coordinate system for the velocity 
vector , we introduce a polar coordinate system (£ p ,0,5 z ) 
which is made nondimensional by scaling with (2RT Q )^ 2 . We 
scale all length dimensions by 'a'. Then the length of the slot 
is 1 i' and its half width 1. We denote by (x,y,z) the scaled 
values of (xl,yl,zl). Thus, 

= 5zl''< 2RT o )1/2 
S 2 = (? xl 2 + 5 yl 2 )/2 RT 0 

and 

9 = Tan ^ (£y-|/£ x i ) 

Finally, we define 

F 1 = i 2 - (2 itRT o )_ 3/2 f (?p> e »5 z ;x,y) . (2.3) 

Then the B6K equation and boundary conditions are as follows: 

S (C0S 0 Sin 6 ly } = 2l<i^ ( )( )1_W (f o _f) ’ 

where 

fo'lMlJ 2 ) 3/2 exp{-B Js.) 
and 

B = [5 z 2 H p 2 +(u 2 +v 2 )/2RT 0 - 2 (u cos 6 + v sin e)y (2RT q ) 1/2 ] 

(2.4) 


wi th , 
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f = exp(-£ p 2 -£ z 2 ) at x = 0, -1 < y < 1 , - J < 9 < \ , 

f = a exp(-£ p 2 -£ z 2 ) at x = 5,, -1 < y < 1 , j <: 0 < , 

? * ■ (2.5) 

v = 0, i.e., /// d£ p d£ z de C p sin 6 f = 0 

and . r ....... 

f = ( f- ) exp(-£ 2 -£ 2 ) at y = 1 , 0 < x < s, 
p o p 

for tt < 0 < 2tt 
and 

f(C p ,0,C z ;x,O) = f(£ p ,-6,£ z ;x,-0) for 0 < x < & . 

The macroscopic quantities are now defined as follows: 

£ = ^ /// % de d ? z C p f 

(pu,pv)/p 0 (2RT 0 ) 1/2 = tt" 3/2 /// d£ p de d£ z (cos e,sin e)£ p 2 f, 

f" = ^’ 3/2 /// d£ de d ? z S D B f > (2.6) 

0 H K 

(P ,P ,P )/p = TT -3 ^ 2 /ff d£ de d£ £ 3 f 
v xx* yy’ xy //K o JJ7 ^p s z s p 

2 2 

• (cos e,sin 6, cos 0 sin e) 

Of great importance is the modified distribution function-- 
g(£p>e;x,y), defined as 

1 00 

g(gp,9;x,y) = -• d- 1 7 2 / d£ f(£ ,e,£ ;x,y) . (2.7) 

Physically, the modified velocity distribution function g(x,y;£p,6) 
is defined to be the number of molecules per unit volume of the 
(x,y) space and the (£ p ,0) space, scaled by the number of 
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molecules per unit volume of the (x,y) space and (£ p ,6) space, 
upstream of the slot. 

The governing equation for g is obtained by integrating 
the B6K equation for f with respect to £ z to give, 


5 P (cos 9 H + sin 9 !y > = m: {r < f" )1 ’“ ( v 9) 


•*-<{& ) 3/2 exp (- Js. 
M o 


The boundary conditions in terms of g are, 

9 = exp(-£ p 2 ) at x = 0, -l<y<l,-|<e<f , 

g = a exp(-£ p 2 ) at x = A, -1 <y<l,|<0<^|, 

v = 0, i.e., // d^ p de £ p 2 sin 6 g = 0 (2.9) 

and 

g = -2- exp(-£ 2 ) at y ■ 1 , 0 < x < a for ir < e < 2ir 
M o v 

and 

g(Sp>e,x,o) = g(^p,-e,x,-o) for o < x < a 


The following macroscopic quantities are definable in terms of 
the moments of g: 
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= // dg_ de E g 
P 0 JJ P P 

(pu,pv)/p 0 (2RT 0 ) 1/2 = // d^ p de ? p 2 g Ceos 6, sin e),-;.0 (2.10) 

(P xx’ P yy’ P xy )/p o = // d? p de V 9 ( cos20 ’ sin20 » cos 0 sin , 0 ) * 

It is evident that the temperature T cannot be defined 
in terms of g. 
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3.0 MOMENT METHOD OF SOLUTION 
3.1 Formulation 

Consider any arbitrary point (x,y) in the slot. A 
discontinuous velocity distribution, made up of four parts— f 1 , f 2 , 
f 3 and f 4 , is assumed at that point and is constituted thus: 

f = f-j for 0 < < oo and 0 < e < tt/2 , 

f = f 2 for 0 < and tt/2 < 6 < it , 

(3.1a) 

f = fo for 0 < E, < oo and tt < 0 < 3tt/2 

o P 

and 

f = f 4 for 0 < C p < 00 and 3ir/2 < 9 < 2u . 

The f . 1 s take the form 
P-; (x »y ) o p 

fi = - exp (-Cp < z ) With i = 1 ,2,3,4 (3.1b) 

The p. 's are unknown functions of (x,y). Any moment 
of f can be calculated in terms of the p^'s. Thus the moment 
<$> of f with respect to a function <J>(£ p ,6,£ z ) is given by, 

oo 00 tt/2 TT 

<4>> = / d? / d? e ( / f-* d0 + / fp<|, de 

-~ o p p o 1 tt/2 ^ 

3tt/2 2tt 

+ / f 3 <|) d 0 + / f.<|> de) . (3.2a) 

tt J 3tt/2 4 


For example. 



(2R ^r/ ' 2 = (pr°2' p 3 +e 4 )/4l,1/2 

The problem is hence resolved if the unknown functions 
p i (i = 1, 2,3,4) are found. These four functions are' determined 
by taking moments of the Maxwell Boltzmann equation. For any 

, , - ! - f w * * & * . • * i ■ . . ■ \ * 1 

quantity «j>(C p » 0 »5 z ) : 

-£ ( HI d '-z % de Sp 2 cos e ! f) L 

+ c /// d 5 z d5 p de c p 2 sin e * f) 

■ ^(r> 1 '“2ra-///<V f >*S d5 P ded5 z 

H 0 0 0 ^ r 

= A g <4>> . (3.3) 

The p. (i = 1 to 4) are replaced by m q 0’ M 01 * M 10 and . M ll as 

unknowns, defined thus: s- 

* 

M 00 = ( p l +p 2 +p 3 tp 4 )/4 

Mq-j = (Pl +p 2"^3" p 4^/^ 7r < 

1/2 ■> ( 3 . 4 ) 

M-| 0 = (p 1 -p 2 -P3 + P 4 )/4tt 7 " 1 

M 11 = ( p 1 “ p 2^ 3 _p 4 ) / 4 TT 

These new variables are noteworthy nsH nee,/- 



18 


p p p 

M = £_ = J. _*£. = 1 = 1 -1L 

00 P 0 2 Po 2 P 0 2 P' 


0 


M 10 - PU/P 0 (2RT 0 ) 1/2 


M„, = pv/p 0 (2RT 0 ) 1/2 


'01 
while 
M 


. W 

n w 0 


The required equations for the four unknowns Mgg, f 

are found by taking <J> = m, m£ cos 9, m£L sin 0 
2 

m^p sin 6 cos 0. The choice of these values for 
by considerations detailed earlier in Section 1.2. 
obtained are as follows: 
continuity equation -- 


9M 


10 


9M, 


'01 


9x 9y 

x-momentum equation -- 


= 0 


1 

2 


9M 


'00 


9M 


11 


9x 9y 
y-momentum equation -- 

9M n i 9M oo 

9x 2 9y 


= 0 


= 0 


while the shear stress equation is 


9M, 


'01 


9M 


10 


9x 

where 


3 y 


(7r) 1/2 M nn T 1-co 
— Tri— * f: > (M lT M 10 -M C 


f - * 1 - ! < M io 2 + M ov >/ M oo 2 


'oi ’ M 10 an ‘ l M 11 
and 

(t> is dictated 
The equations 


(3.5) 


/ M 00^ ’ 
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The boundary conditions expressed earlier in terms of 
f can be expressed now in terms of ^go’ M 01’ ^10 anc * ^11* They 
are as follows: 


At x = 0, 

+ tt 1/2 M, 


and 


'00 " "10 


" 1/2 H 11 + M 01 


= 0 


At x = a. 


M - tt ^ 2 M 
n 00 1 0 


and 


„ 1/2 M n - M 0) 


= 0 


At y = 1 and y = 0, 


M 01 - 0 


and 


* 1/2 M n - M 10 “ 0 


(3.6) 


The above system of equations is hyperbolic. The boundary 
conditions are, however, unorthodox since only two (p-| and p. ), of 
the four unknown quantities are specified at x = 0; the other 
two (pg and p^) are specified at x = l. This is a boundary 
value problem, and to ^reduce it to a classical initial value prob- 
lem, all four quantities p. have to be specified at an initial 
front, say x = 0. However, the question arises: Does a unique 

solution exist for the above system of equations with those 
unorthodox boundary conditions? Using the work of Saranson (1962), 
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it is verified that a unique solution does exist for all values of 
the parameters. The details are elaborated in Appendix A. 

3.2 Free Molecule Solution 

i ' 

In the limit Kn Q -*■ °°, which corresponds to the free 
molecule flow, the above equations can be solved in terms of two 
unknown functions G and H, giving 

M oo = G ( x+ y) + G ( x -y) 

M n = - \ (G(x+y) - G(x-y)> 

(3.7) 

M oi = H ( x -y) ■ H (x + y) 

M 1q = H(x+y) + H(x-y) 


In general the corners of the boundary are points of 
singularity. If we assume that at the corners of the boundary 
(x = 0 and y = 0 , y = 1 ; x = £ and y = 0 , y = 1 ) multiple 
boundary conditions can be simultaneously applied, then closed 
form solutions can be obtained for the G and H functions for 
any rational value of £. However, the expressions for G and 
H for any rational value of £ are considerably complicated. 
Only integer values of £ are considered and the results for G 
and H are tabulated in Appendix B. An interesting feature of 
these results is that the solutions differ depending on whether 
£ is an even integer or an odd integer. Of particular interest 
is M 1q , which is the mass flux pu scaled by ( 2RT Q ) 1 p Q . 
When £ is an even integer, 



. ( 1 -«? 

2Tr 1/2 [l+(£/7r)] 


(3,8) 
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When Z is an odd integer. 


M 


10 


0-g) 


o 1/2/, . & 1 \ 

2,1 (1 + ? • 7tI+2T 1 


(3.9) 


Both the above expressions are valid for all values of x and y 
in the slot. Since M-jq is a constant everywhere in the slot, 
hence 2M 1Q gives the total mass flux through the slot scaled by 

P 0 < 2 RT o > 1/2 a - 

A discrepancy in the two expressions is evident. For 
large lengths, it is negligible, but becomes significant for shorter 
lengths. In Fig. 4 a plot of the mass flux through the slot vs. Z 
has been made and is compared with the accurate numerical solution 
of Reynolds and Richey. Our solution for the mass flux is lower 
compared to Reynolds and Richey's solution. This error, for the 
range 1 to 12 for £, increases as Z increases, being as much 
as 35% for Z - 12- In fact, if our calculations were continued 
for Z greater than 12, the .errors should get progressively worse. 
Also the difference in expressions for even and odd lengths is 
brought out as ripples in the curve (see Fig. 4). This ripple, 
however, dies out as the length Z increases beyond Z = 6. The 
inevitable conclusion is that the solution is significant only for 
Z between 6 and 12; beyond Z = 12, errors involved in the calcu- 
lated value of the mass flux get very large as Z increases. 

. ; ' * 1 " 

3.3 Transition Flow Calculations 

The system of hyperbolic equations is solved by the 
method of characteristics. 
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The characteristics are given by: 



(3.10) 


The compatabil ity conditions along the characteristics 


are: 


dy. 

dx 


d£ 

dx 


dM, 


'00 


dM 


dx 


+ 2 


11 


dx 


= 0 


dM 


10 


dM 


01 


TT 


1/2 


(3.11) 


dx 


dx 


Kn 


owWio> ■ 0 


The grid points are taken as shown in Fig. 2. NN grid 
points are taken at the entrance, spaced at a distance 1 /NN from 
each other. The grid points immediately adjacent to the boundaries 
(the centerline and the wall) are at a distance 1/2NN from the 
boundaries. The two characteristics dy/dx = ±1 pass from each 
grid point, forming a grid pattern, as shown in Fig. 2. The 
length i is always so chosen that 2i,NN is an even integer. 

Integrating the above equations and expressing the 
results in a finite difference form, we get the following expression 
for any three points A, B and C (see Fig. 2): 

Along &-1. 

M oo (A) + 2 M 11 (A) = M oo( B) + 2 M 11 (B) 

and (3.12) 

M 10 (A) + M 01 (A) = M 10 (B) + M 01 (B) 

1/2 

- <« M ooWio>< t; 

+ f< M oo M ir M oi M io>< r )’■“>) 

o B 
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Alcg 

M oo (A) ■ 2 M n (A) = M oo (c) ■ 2 M n (c) 

and 


- M 1q (A) + M q1 (A) = - M 1q (C) + M q1 (C) 


(3.13) 


with 


T 

T7 


= 1 - 


1/2 

%r— (t(M oo M n' M oi M io )( {" > > 


o C 

2 ( M 10 2+M 01 2 > 

3 


M, 


00 


and 


Ax = 1/2NN 


Thus knowing the various M's at B and C, . the M's 
at A can be immediately obtained. Thus if all the M's are 
known at x = 0, then the M's at all other grid points can be 
obtained with ease. When w ^ 1 , at each point an iteration is 
involved to calculate M-|q and Mqi . When u) = ,l, an iteration 
is involved only to calculate M-jq at the wall. Define 


P' 

Q' 


M oo - M 10 

J' 2 M„ - M 01 


R' 

S' 


M oo + ul/2 M 10 


* 1/2 M 11 + M 01 


(3.14) 
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A survey of the boundary conditions at x = 0 reveals 
that R' = 1, S' = 0. No values for P' and Q 1 are available 
at x = 0 to proceed with the initial value computation detailed 
above (P ' and Q' are available only at x = £, where P' = a, 

Q' =0). However, an iteration procedure detailed in Appendix C 
can be followed to arrive at the right starting values which would 
give P' = a, Q‘ = 0 at x = £. For the range of a, Kn Q and 
£ considered, 12 was the maximum number of such iterations required. 

A Fortran program was written to solve the system of alge- 
braic equations. The following are the parameters in the problem: 
w, £, NN, Kn Q and a. For the range of parameters considered, it 
was found that T/T q differed from unity only by 0(10~ 2 ). 

Because of this weak influence of T/T Q , it was decided to keep 
w = 1 and eliminate the iteration, except at the wall. For Kn Q > 1, 
NN could be kept at a nominal value— between 3 and 5. However, when 
Kn 0 < 1» if NN was kept too small, the iteration scheme for M-jq 
at the wall did not converge. The following criteria were evolved 
for NN (for Kn Q < 1) so that the iteration scheme for M-|g at 
the wall would be convergent: 

NN Kn Q > (tt)~ 1/2 - 0.25 a . (3.15) 

Kn Q was given the values 5, 1 and 0.5, while £ had values 1, 

2, 4, 8 and 12. 

The practice evolved in the computation was to keep £ 
and a fixed and vary Kn Q — starting with Kn Q = 5 and decreasing 
it to Kn Q =0.5. To hasten the initial value iteration process. 
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values obtained at x = 0 for P 1 and Q' at the previous Kn Q 

were used for the next Kn . 

o 
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4.0 DISCRETE VELOCITY ORDINATE METHOD OF SOLUTION 

The discrete velocity ordinate method of solution takes cogni- 
zance of the fact that the distribution function (modified or 
otherwise) is a function of both space and velocity variables. 
However, in this method only certain discrete values of the 
velocity are considered and the distribution function (modified or 
otherwise) is evaluated at these discrete values of the velocity. 
4.1 Formulation 

In this method, the modified distribution function 
g(£p»6;x,y) is the fundamental quantity that is evaluated. The 
function g, defined by (2.7), is governed by Eq. (2.6) with the 
boundary conditions (2.10). 

Once g is known, most moments, except for T/T , can 
be immediately evaluated— see Eqs. (2.10). A survey of the moment 
solution shows that (T/T ) varies from unity by less than 0.5% 
for the slot lengths considered. In consequence (T/T ) is set 
identically equal to unity. Thus, Eq. (2.8) has the simplified 
form, 

y cos e|f + sin e |l) =^-^-(90-3) (4.1) 

where 

g 0 = ^-exp{-(£ p cos 6 ■ u ') 2 - (£ p sin © - V') 2 } , 

u 1 = u/(2RT p ) 1/2 and v' = v/(2RT Q ) 1/2 

The boundary conditions (2.9) remain the same. 

Formally, g is a function of four independent, con- 
tinuous variables, x, y, 6 and £ . However, the governing 
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equation for g, (4.1), has no derivatives with respect to either 
£p or 0. Hence, in accord with the philosophy of the discrete 
ordinate method, g is considered only at certain discrete values 
of £p and 6. The choice of these values is strongly motivated 
by the physics of the problem. 

Considering the present method as a natural extension of 
the moment method, the eight discrete values of 9 are given by 
9 = with M = 0 to 7. The choice of the discrete values of 
£p is dictated by the consideration that our final interest is not 
in g itself but in the moments of g--most notably the mass flux 
pu along the slot. For a fixed 9, the free molecular and 
continuum functional form for g is, 

exp{-(£p cos 9 - u') 2 - (5 p sin 0 - v') 2 } . 

It seems reasonable to choose a quadrature formula that will be 
accurate for this form of g. As u' and v 1 are usually « 1, 

we linearize accordingly. The evaluation of p, pu, pv and P. . 

(i,j = x,y) involves integrals, with respect to E, , of the form 

«. -£ 2 
{ p VV 

where h^p) is a polynomial of degree 3 or less. We use a 
Gaussian quadrature formula for the evaluation of such integrals. 
Then the two discrete values of £ chosen are the Gaussian roots 
corresponding to the above integrals, denoted by £ a and £ b . Let 

the corresponding Gaussian weights be denoted by w a and w^. The 

numerical values of E, & , £ b , w a and w b are evaluated in Appendix 
D. 
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The moments of g, besides being the quantities of primary 
interest, are, as seen from Eq. (4.1), involved in the solution of 
g itself. These moments involve integration of g with respect to 
both £p and 6. As mentioned earlier, the integration with 
respect to 5 is done using the Gaussian quadrature formulae. 

This is done for all values of Kn Q , and it is anticipated that 
the accuracy of such a scheme would increase as Kn Q 0 and °° 

p 

when the speed distribution will be proportional to exp(-£p ). As 
for the integration of g with respect to 0, there are some 
important aspects to the flow that have to be considered. When the 
flow everywhere is close to continuum, g varies continuously with 
respect to 0. As the flow tends further and further from continuum, 
g varies more and more steeply with respect to 0. Finally, in the 
limit of free molecule flow, g varies discontinuously with 
respect to 0. To see this, consider an arbitrary point (x,y) in 
the slot. Then the distribution function g in the solid angles 

p 

subtended by the entrance and the exit are exp(-£p ) and 
2 

a exp(-£p ), respectively. The number flux from the wall varies 
almost linearly along the length of the slot (see Reynolds and 
Richey, 1967). Translating this, the distribution function of the 
particles coming from the wall is close to the form (a(x,y) + 

p 

cot 0 b(x,y)}exp(-£p ), where 'a' and 1 b 1 are functions of 
(x,y) and 0 lies in the solid angle subtended by the wall at 
(x,y). 

First, a simple trapezoidal rule is used for integrating 
over the 0 space, using the eight discrete values in the 0 space. 
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This scheme is very good when the g varies slowly with respect 
to 0. However, as the flow tends to the free molecule regime, 
g varies steeply and finally discontinuously with respect to 6. 
Here this method of integration is inadequate. 

The second mode of 6 integration allows for the dis- 
continuous form of g (with respect to 0). In the solid angle 
subtended by both the entrance and the exit, the trapezoidal 
integration is used. In each of these solid angles, 1 to 3 dis- 
crete ordinates can lie in their span for each of the discrete 
values of £ p . In the solid angle subtended by each of the walls, 
the integration is done assuming that g varies with 0 as 
cot 0. In each of the solid angles subtended by the wall, for 
each discrete ordinate of £ p , 2 to 3 discrete ordinates in 0 

can lie in their span. Thus, if there are two such ordinates we 
assume that 

g(x,y;Cp,e) = {a(x,y;£ p ) + b(x,y;£ p )cot e} exp(-£ p 2 ) (4.2a) 

and evaluate the two unknown functions a(x,y;£ p ) and b(x,y;£ p ) 
from the known values of g for the two 0 ordinates, for each 
speed £ p . If, however, 3 discrete ordinates lie in the solid 
angle, assume that 

g(x,y;£p,e) = (a(x,y;£ p ) + b(x,y;£ p )cot 0 + C(x,y;£ p ) 

(|i-e) 2 exp(-s p 2 ) , (4.2b) 

where k = 3,1 for the top wall and bottom wall, respectively. 

The unknown functions a,b,c are evaluated from the known values 
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of g at the three e ordinate values (for each speed £p). 

This method of integration is done not only for the free molecule 
flow, but also continued into the transition regime. 

As emphasized earlier, the trapezoidal integration is 
increasingly accurate as the flow tends away from the free molecule 
regime. The cot 0 integration is certainly accurate in the free 
molecule regime; for it to be meaningful in the transition regime, 
the results obtained by this method have to tend to the results 
obtained by the trapezoidal integration. 

Let Q(£ ,0) = Q, (£ ) Q 2 ( 0 ) * 3e a separable function 
whose moment (with respect to g) is desired. Thus, by the 
trapezoidal integration in 0, 

CO 2tt 

/ / d£ d0 £ Q(£ ,0) g(x,y;£ ,0) 

0 0^^^ H 


1 1 W. Q-i(?-j) exp(£- ) g(x,y;£. , —r ) Q 2 ( -jr ) 

i=a n=0 


(4.3a) 


By the cot 0 integration scheme for 0, 

co 2ir 

/ / d£ d0 £ Q(£ ,0) g(x,y;£ ,0) 

0 0 ^ ^ H " 


= J w. Q^S-j) exp(£ n . 2 ) / Q 2 (0) g(x,y;£ ,0) d0 

i=a o ^ 

b 2 e l 

= I w i (€4 ) exp(£. ) { / d0 Q ? (0) g(x,y;£, ,0) 
i = a 1 e c 1 

e 2 

+ / d0 Q 2 (e)[a 1 (x,y;£ i ) + b.|(x,y;£.) cot 0 + 

9 1 
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+ c 1 (x,y;c i )( \ - e) 2 ] + / de Q 2 (e) gfx.yi^'.e) 

9 2 

9 4 • • . 

+ / de Q 2 ( 0) [^2^ x ; ^i ) + b 2( x >y ; ^i^ cot 6 

e 3 ' 

+ c 2 (x,y;^ i )( - 9) 2 ] > > (4.3b) 

where 

e 0 - tan’ 1 ( If ) . 

6 1 = Tan -1 ( ) 

©2 = tt - Tan 1 ( ) 

e 3 = it + Tan" 1 ( ) 

6 4 = 2it - e o • * 

The first and third terms in the { } bracket represent integra- 
tion over the slot entrance and exit solid angles, and are 
evaluated by trapezoidal integration; the second and fourth terms 
represent integration over the solid angle of the bottom and top 
wall. Substituting Q(5 ,0) = 1, 5 p cos 9, £ p sin 9, etc., we can 
easily evaluate p, pu, pv, etc., from both (4.3a) and (4.3b). 

The grid points are taken as shown in Fig. 3. NN grid 
points, equally spaced, are taken in the y direction starting 
from the wall and ending at the centerline. The index J is used 
for the grid points in the y direction'; with J =1 representing 
the wall grid point and J = NN the centerline of the slot. In 
the x direction, the grid points are equally spaced by a distance 
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1/(NN-1); they start at the entrance and end at the exit. The 
index I is used for the grid points in the x direction with 
1 =1, representing the grid points at the slot entrance and 
I = N = (NN-1) &+1 , the grid points at the slot exit. The total 
number of grid points to be considered equals NN*N. The grid 
points at the slot entrance and exit could be considered either 
just inside or outside the slot. The boundary conditions at the 
slot entrance and exit are accordingly different. 

Let us use the notation g(I,J,K,n) to denote the value 

of g at the grid point located at (x,y) given by 

x = (1-1 )/ ( NN-1 ) 
y = (J-l )/ (NN-1 ) 

K = 1 and 2 are the respective values used to denote 5=5 

P a 

and 5p = 5^* while n can have the values 0 to 7, and this is 
defined by 0 = mr/4. Consider the case where the grid points at 
the slot exit and entrance are just outside the slot. Then the 
boundary conditions at the slot entrance and exit grid points are 
as follows: At the slot entrance g(l,J;K,n) = exp(- 5 p 2 ), for 

J = 1 to NN, K = 1 and 2 and n = 0,1 ,2,6 and 7. In addition, 

at the slot corner (x = 0 , y = 1 ), it is assumed that 
9(1 .1 ;K,5) = exp(-5p^) for K = 1 and 2. At the slot exit, 
g(N,J;K,n) = a exp(-5 p 2 ) for J = 1 to NN, K = 1 and 2, and 
n = 2 to 6 . Also at the slot corner (x = 1, y = 1) it is assumed 
that g(N ,1 ;K,7 ) = a exp(-5p 2 ) for K = 1 and 2. 

Consider the case where the grid points at the slot 
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entrance and exit are taken just inside the slot. Then the boundary 
conditions at the slot entrance and exit grid points are as follows: 

At the slot entrance g(l,J;K,n) = expC-f^) for J = 2 to NN, 

K = 1 and 2, and n = 0,1 and 7. At the wall, J = 1, g(1,I;K»n) 

= exp(-£p ) for K = 1 and 2 and n = 0 and 7; also the reflec- 
tion from the wall is diffuse and v = 0. At the centerline, J = NN, 
g ( 1 ,NN;K,n) = g(l ,NN;K,8-n) for K = 1 and 2, n = 1,2,3. At the 
slot exit g(N,J;K,n) = a exp(-£p^) for J = 2 to NN, K =1 and 2, 
and n = 3, 4 and 5. At the wall, J = 1, g(N,l;K,n) = a exp(-£p ) 
for K = 1 and 2, and n = 3 and 4; also the reflection from the wall 
is diffuse and v = 0. At the centerline, J = NN, g(N,NN;K,n) = 
g(N,NN;K,8-n) for K = 1 and 2, n = 1, 2 and 3. 

In addition to the above two possible boundary conditions 
at the slot entrance and exit we have the following boundary condi- 
tions valid for I = 2 to (N-l). At the wall, J = 1, v = 0 and 
the reflection from the wall is diffuse. At the centerline, J = NN, 
by symmetry about the y axis, g(I,NN;K,n) = g(I ,NN;K,8-n) for K = 1 
and 2 , and n = 1 , 2 and 3. 

Numerical experiments with the above two possible sets of 
boundary conditions (at the slot entrance and exit grid points) were 
performed. The free molecule results obtained by assuming that the 
grid points are just outside the slot entrance and exit are much 
better (by comparison with the solutions of Reynolds and Richey). 
Henceforth it is assumed that the grid points at the slot entrance 
and exit are just outside the slot. 
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Checking our list, we find that with the grid points at 
slot entrance and exit placed outside the slot, (6NN-4) quanti- 
ties have to be assumed to specify g(l,J;K,n) completely for 
all values of J, K and n. In the case of free molecule flow 
(see Section 4.2) several simplifying features exist. Conse- 
quently, there is a rather dramatic decrease in the number of 
quantities to be assumed to specify g(l,J;K,n) completely for 
all J, K and n. 

At the wall we note that g(I,l;K,0) and g( I ,1 ;K,5) 
are both discontinuous and have two values each. This is because 
of the two faced character of the distribution function g at 
the wall, since there are essentially two totally different 
classes of molecules there— the molecules reflected from the 
wall and those coming into the wall. 

4.2 Free Molecule Limit Solution 

In the free molecule limit, i.e., Kn Q <», Eq. (4.1) 
reduces to, 

cos 6 |f + sin 0 = 0 (4.4) 

Integrating (4.4) along the characteristic direction, dy/dx = 
tan 6, the solution of (4.4) is that g is a constant along the 
characteristic direction dy/dx = tan 0, for all £ p . 

In terms of our grid points set and discrete velocity 
ordinates this can be interpreted as. 
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g(I+l,J;K,0) = g(I,J;K,0) 
g(I+l,J;K,l) = g(I ,J+1 ;K,1 ) 
g(I+l,J;K,2) = g(I+l ,J-1;K,2) 
g(I+l,J;K,3) = g(I,J-l;K,3) 

, (4.4a) 

g(I+l ,J;K,4) = g(I,J;K,4) 

g(I+l ,J ;K,5) = g(I,J+l;K,5) 
g(I+l ,J:K,6) = g(I+l ,J+1 ;K,6) 

g(I+l ,J;K,7) = g(I ,J-1 ;K,7) 

for 1=1 to N-l, K = 1 and 2 and J = 1 to NN. 

To reduce the two point boundary value problem to an 
initial value problem, it is required that g(l,J;K,n) be com- 
pletely known for J = 1 to NN, K = 1 and 2 and n = 0 to 7. 

As discussed in an earlier section, in general (6NN-4) of these 
quantities have to be assumed. However, for this free molecule 

flow situation the distribution function is, for 0 constant, 

o 

directly proportional to exp(-£p ). Hence only K = 1 has to 
be considered. Also, for 0 = tt , i.e., n = 4, the value is 

given by the value at xl = Jla, i.e., I = N. Hence we need to 
assume only (2NN-2) quantities at xl = 0, i.e., 1=1. 

When g(l,J;K,n) has been assigned for J = 1 to NN, 
n = 0 to 7 and K = 1, say, the values of g(2,J;K,n) can be 
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computed as follows: 

At the wall we have from (4.4a), 

g(2,l ;K,5) = g(l,2;K,5) 

i.e., we know the distribution function for molecules leaving the 
wall at 0 = 5 tt/ 4. By the condition of diffuse reflection, 

g(2,l;K,n) = g(2,l;K,5) for n = 4, 6, 7 and 0, 

with the understanding that the values for n = 0 and 4 represent 
the outgoing distribution function. We now have g(2,J;K,6) = 
g(2,J;K,2) for J = 1 to NN by symmetry (the number flux at the 
two walls is the same for the same x value). Equations (4.4a) 
now serve to determine g(2,J;K,n) for all values of J = 2 to 
(NN-1), and g(2,NN;K,n) for n = 4, 5, 6 and 7, i.e., molecules 

crossing the centerline from y > 0. By the symmetry at the center- 
line we also have, g(2,NN;K,l ) = g(2,NN;K, 7) , i.e., 0 = ir/4 and 
7 tt/ 4 and g(2,NN;K,3) = g(2,NN;K,5). The only remaining value is 
that for the wall and 0 = 3 tt/ 4, i.e., g(2,l;K,3). While all 

previous relations are exact for free molecule flow we have to 
determine g(2,l;K,3) by the condition of zero normal mass flux 
at the wall. The accuracy of this calculation depends on the accu- 
racy of the 0 quadrature formula. 

The above process can be repeated for all values of I 
up to N. At I = N the calculated values of g(I,J;K,n) for 
J = 1 to NN and n = 3 and 5 have to be equal to a exp(-Cp^). 

If they are not, an iteration scheme, similar to that detailed in 
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Appendix C, can be used to obtain those values of g(l,J;K,n) 
for J = 1 to NN and n = 3 and 4. which give g(N,J;K,n) = 
a exp(-£ *) for J = 1 to NN and n = 3 and 5. A total of 

r 

(2NN-1) sweeps from 1=1 to N are required for each such 
iteration. 

The solutions for g(I,J;5p,0) were obtained by both 
the trapezoidal and cot 6 integration schemes. These solutions 
are different as the computational form for the boundary condition 
pv = 0 at the wall is different by the two methods. The macro- 
scopic quantities and the mass flux through the slot (by trape- 
zoidal integration of pu across the slot section) were also 
obtained by both forms of 0 integration. The only parameter 
varied is the slot length i,, since the free molecule solution 
depends linearly on a. 

I 

The mass flux through the slot obtained by the trape- 
zoidal integration (with respect to 0) agrees well with the 
results of Reynolds and Richey for small slot lengths. However, 
as the slot length increases there is a systematic overestimation 
of the mass flux with an error of 25 % for a slot length of 12. 

By the cot 0 integration, there is a spectacular agreement with 
the mass flux results of Reynolds and Richey with an error of less 
than 1% for a slot length up to 12. This is not surprising as we 
took advantage of the linear wall flux dependence to pick our 0 
quadrature scheme. 

Figure 4 displays the mass flux through the slot, as a 
function of slot length, by both the cot 0 and trapezoidal 
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integrations. The slot lengths were considered over the range 
1 to 12. 

4.3 Transition Flow Calculations 

Here the entire equation (4.1), with associated boundary 
conditions for g has to be considered. Replacing the left hand 
side of Eq. (4.1) by the derivative along the characteristic direc- 
tion dy/dx = tan 0, dividing by and introducing a quantity L, 

a? • L ■ £n t aa : (9 °~ 9) • (4 - 5) 

p K 0 0 

where dx/ds = cos 0 and dy/ds = sin 0. Expressing (4.5) in a 
simple finite difference form, we have 

g(s+As;? p ,0) = g(s;£ p ,0) + L(s,? ,0) As , (4.6) 

where L, the right hand side of Eq. (4.5), is evaluated at s. 
Equation (4.6) can be specialized for the discrete velocity ordi- 
nates values. Let us again use the notation g(I,J;K,n) to denote 
the value of g at the grid point (I,J) with K = 1 and 2 
representing £ p = and £ b , while 0 = nir/4. Then, for our 
grid points, the following equations are valid: 

g(I+l,J;K,0) = g(I,J;K,0) + L(I ,J;K,0) 

9 l/2 

g(I+l,J;K,l) = g(I,J+l;K,l) + L(I,J+1;K,1) 

g(I+l ,J;K,2) = g(I+l,J-l;K,2) - ppyy L( 1+1 ,J-1 ;K,2) 

2 l/2 

g(I+l,d;K,3) = g(I,J-l;K,3) - L(I ,J-1 ;K,3) (4.7) 

g(I+l,J;K,4) = g(I,J;K,4) - XNN^TT l (I,J;K,4) 
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, 1/2 

g(I+l,J;K,5) = g(I,J+l;K,5) - -pppy L(I,J+1;K,5) 

g(I+l ,J;K,6) = g(I+l ,J-1 ;K,6) + L(I+1 ,J-1 ;K,6) 

g(I+l,J;K,7) =g(I,J-l;K,7) +-^ rr yL(I,J-l;K,7) 

where I = 1 to (N-l), J = 1 to NN and K = 1 and 2. 

For purposes of computation, g(I,J;K,7) has to be 
completely specified at 1=1 for J = 1 to NN, n = 0 to 7 
and K = 1 and 2. As detailed in Section (4.1), (6NN-4) 

quantities have to be assumed to fully specify g(I,J;K,n) at 
1=1 for all values of J, K and n. 

For the transition flow, there is no simple relation con- 
necting g ( I ,J ;K,n) for different values of K. Once again, for 
I = 2 to (N-l) and K = 1 and 2, both g(I,l;K,0) and g(I,l;K,4), 

1. e., the value of the distribution function at the wall for 

6 = 0 and tt is double valued. Also, we find that for K = 1 and 

2, g(l,J;K,4), i.e., 9 = ir is not immediately known. Further, 
we find that g(I,l;K,3) for K = 1 and 1=2 to (N-l) have to 
be assumed. Thus in all (6NN+N-6) quantities have to be assumed 
in all to reduce the two point boundary value problem to an initial 
value problem. Therefore, the transition flow calculations are 
considerably more complicated than the free molecule flow calcu- 

1 ati ons . 

When g(l,J;K,n) has been assigned for J = 1 to NN, 

K = 1 and 2 and n = 0 to 7, the values of g(2,J;K,n) can be 
computed as follows: 
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At the wall, using (4.7), we can calculate g(2,l;K,5) 
in terms of g(l,2;K,5) for both K = 1 and 2, i.e., we know the 

distribution function for molecules leaving the wall at 0 = 5ir/4. 

By the condition of diffuse reflection, 

g(2,l;K,n) = g(2,l;K,5) for n = 0,4,6 and 7 and K = 1 and 2. 

The understanding here is that the values for n = 0 and 4 repre- 
sent the outgoing distribution function. Further, using Eqs. (4.7), 
g(2,l ;K,n) can be calculated for K = 1 and 2 and n = 0,1 and 4.; 
the values here for n = 0 and 4 are understood to represent the 
incoming distribution function. Since g (2 ,1 ;1 ,3) , i.e., the 

value of g at the wall for E = E and 9 = 3 tt/ 4, has been 

P « 

assumed, hence all values of g(2,l;K,n) have been found except 
for g(2,l;2,3) and g(2,l;K,2) for K = 1 and 2. These quanti- 
ties correspond to £ p = C b , 8 = 3tr/4 and £ p = £ a ,S b> e - tt/2, 
respectively. As an initial guess for g(2,l;K,2), assume that 
g(2,l;K,2) = g(2,l;K,6), i.e., g is the same for 9 = tt/2 and 
9 = 3tt/2, for K = 1 and 2. Then g(2,l;2,3) can be calculated 
using the boundary condition pv = 0 at the wall. Thus g(l,J;K,n) 
is completely specified at the wall, J = 1. 

Using (4.7), g(2,J;K,n) can now be completely specified 
for J = 2 to (NN-1), K = 1 and 2 and n = 0 to 7. At the center- 
line J = NN, g(2,NN;K,n) can be calculated for K = 1 and 2 and 
n = 0,2, 3, 4, 6 and 7. Symmetry conditions there also give. 


g ( 2 ,NN ; K ,5 ) = g(2,NN;K,3) 
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and 

g(2,NN;K,l ) = g(2,NN;K,7) for K = 1 and 2. 

Further, symmetry conditions at J = NN require that g(2,NN;K,2) = 
g(2,NN;K,4) for K = 1 and 2. If our calculations do not give 
this, it means that the assumed values of g at the wall for 
9 = tt/2 and E = E and £. , i.e., g(2,l;K,2) are not appropriate. 

pa D 

An iteration scheme similar to that in Appendix C can be followed 
to obtain accurate values of g(2,l;K,2) for K= 1 and 2 which 
would give, to sufficient numerical accuracy, g(2,NN;K,2) = 
g(2,NN;K,6). Thus, g(2,J;K,n) can be completely calculated for 
all values of J, K and n. 

The above process of calculation used for g(2,J;K,n) 
can be repeated for all values of I up to N. The check to see 
if the (6NN-6+N) values that were assumed were the appropriate 
ones consists of two parts. First, since the reflection from the 
wall was assumed to be diffuse, it is required that 

g(1 ,1 ;1 ,5)exp(£ g 2 ) = g(1 ,1 ;2 ,5)exp(£ b 2 ) for I = 2,(N-1). 

This gives (N-2) matching conditions. Further, it is required 
that g(N,J;K,n) = a exp(-£p^) for J = 1, NN, K = 1 and 2 and 
n = 3, 4 and 5. If these conditions are not satisfied an iteration 
scheme like that in Appendix C can be used to get the (6NN-4+N-2) 
values required. For each such iteration (6NN+N-5) sweeps are 
required through the field 1=1 to I = N. 

Calculations were performed for % = 1,2,4,8,12, a = 0.8, 

0. 5,0.1 and Kn Q = 5. 0,1. 0,0. 5. The procedure used was to keep 
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the length % constant first. Then a is systematically varied. 

However, for each a, Kn Q is systematically varied from 5.0 to 

0.5. The initial guess given for Kn Q = 5.0 is taken from the 

free molecule flow values. The initial guess for subsequent values 

of Kn„ is the calculated value at the previous Kn . The 0 
o o 

integrations were performed by both the trapezoidal and cot 0 
integrations. For the parameters considered (Kn Q < 5) the agree- 
ment by the two modes of 0 integration is very good, though the 
computer time involved for the cot 0 integration is considerably 
more than for the trapezoidal integration. Also as the slot length 
is increased, the computer time consumed in the calculation in- 
creases as the square of the slot length. For a slot length of 
12 a number of minutes of computer time (CDC 6400) are required. 
This could be pinpointed as one of the deficiencies of the method. 

A comparison of the solutions obtained here with the moment method 
shows a remarkably good agreement— the difference is less than 5 % 
for Kn Q = 0.5, 1.0 and 5.0. However, as Kn Q gets larger, a 
considerable difference in solutions is expected between the two 
methods. This is because the free molecule solutions by the two 
methods are so different. 
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5.0 DISCUSSION OF SOLUTIONS 

Solutions by both the moment and the discrete ordinate methods 
have been obtained for free molecule and transition flows. The 
latter method employed two forms of angular integration in the 
velocity space, the trapezoidal rule and the cotangent integration. 

The range of parameters considered is a = 0.1, 0.5 and 0.8, 
l = 1, 2, 4, 8 and 12, and Kn Q = <», 5, 1.0 and 0.5. 

5.1 Free Molecule Results 

The free molecule solutions are compared with the solu- 
tions of Reynolds and Richey (1967). The number flux from the wall 
by the moment method is shown in Figs. 7 and 8 for lengths l = 4 
and 8; Figs. 5 and 6 show the same by the discrete ordinate method 
with both forms of integration. Surveying these figures, we find 
that the discrete ordinate method with the cotangent integration 
agrees closely with the linear wall flux variation obtained by the 
solution of Reynolds and Richey. This is not surprising, since we 
had specifically tailored our 6 integration to a linear wall flux 
variation. The wall flux by the moment method and the discrete 
ordinate method with a trapezoidal integration seem to follow the 
overall linear variation of Reynolds and Richey but have step like 
variations. These steps correspond to sudden changes in the wall 
flux values which occur for even integer values of x, the coordinate 
along the slot axis. This is attributable to the fact that the only 
oblique angular direction considered in the velocity space is ir/4 
with respect to the slot axis, and hence there is a certain periodi- 
city over distances equal to the total slot width (= 2). Figure 4 
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shows the mass flux through the slot by the various methods. The 
moment method underestimates the mass flux; the error continues to 
increase as the slot length increases. On the other hand, the 
discrete ordinate method with the trapezoidal integration over- 
estimates the mass flux; the error continues to increase as the 
slot length increases. In both these cases, the poor results are 
due to no allowance being made for the discontinuity of the distri- 
bution function with respect to 6 and also for the linear varia- 
tion of the distribution function along the slot axis. The discrete 
ordinate method with a cotangent integration makes allowance for 
both these factors. Hence the mass flux computed here agrees very 
well with that by Reynolds and Richey. 

5.2 Transition Flow Results 

In the transition regime, for the values of Kn Q con- 
sidered, the distribution function varies continuously with respect 
to 6 (except at the walls). Hence the differences in the discrete 
ordinate solutions with trapezoidal and cotangent integrations is 
barely discernable. This can be seen, for example, in Fig. 9, where 
the wall flux values are plotted by the two methods for l = 1 , 
a = 0.1 and Kn Q = 0.5. These differences by the two methods should 
get less as the slot length increases. For the range of parameters 
considered the difference in the computed mass flux through the slot 
is less than 4%. Figure 10 shows a plot of the number flux from the 
wall for i = 4, a = 0.1, Kn Q = 5.0 and 1.0 obtained by using the 
discrete ordinate method with trapezoidal integration. The number 
flux seems to vary slowly with respect to x except close to even 
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integer values of x. The significance of the even integer values 
of x seems the same as for free molecule flow. 

The moment method calculation in the transition regime 
yields results which compare very favorably with the discrete 
ordinate method solution. Figure 11 shows the plot of the number 
flux at the wall. The pattern seems to match that of the discrete 
ordinate solution, although the variations at even values of x 
seem sharper. The significance of the even values of x is the 
same as in the free molecule flow. Figure 12 shows the number flux 
from the wall vs. x, by the moment method for l = 12, a = 0.1 
and Kn Q = 5.0, 1.0 and 0.5. Remarkably, any periodicity with 
respect to even values of x is absent. Further, the variation 
is almost linear. Also the values of the mass flux obtained by the 
moment method agree very well with that by the discrete ordinate 
method. For the range of parameters considered, the difference 
seems to increase as Kn Q increases, being a maximum of 2% for 
Kn Q =0.5 and a maximum of 6% for Kn Q = 5.0. Further, the plane 
Poiseuille formula with slip boundary conditions (valid strictly 
only for an °° length slot) gives values for the mass flux through 
the slot that agree with the moment and discrete ordinate solutions 
for Kn Q = 0.5 and 1.0. The errors are less than 5 %. The remark- 
able feature here is that this result is valid for all lengths 
considered Z = 1, 2, 4, 8 and 12. However, we have to substitute 
the actual values of the average pressure at the slot entrance and 
exit rather than the values p Q and ap Q , respectively. 

For the parameters considered, the following interpolation 
formula is proposed for m, the mass flux through the slot: 
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m 


_ 1 , 0.1331 


Kn. 


= 


Kn. 


p 0 a(2 RT 0 ) 1/2 1 -o 

Kn. 


0+a)T-| + 2(T-| - -gS-log S ] ) 


+ — z lo 9 h iog( - * j>* ~ r ‘ t ) > 

8/F 1 2 


v4, 2 +4 


where 


T - (1 "0Q& 

T 1 ” A + 4.85 - 2.85 exp{0.21 (Kn Q -5)} 


S, = 


2Kn„ + 8 ( 1 +a+T ) 


1 2Kn 0 + 8(l+a-T) 


The numerical values of m obtained by calculations differ from 
this formula by less than 5%. The above formula is patterned after 
the formula proposed by Fryer (1966) for the mass flux through an 
infinite cylindrical tube in the transition regime. Crudely, in 
the expression for m, the first term is the mass flux in the 
continuum regime, the second term the contribution due to slip at 
the boundary, and the last term the free molecule mass flux. 

A systematic study on the effect of the number of grid 
points NN on the results obtained, both by the moment and discrete 
ordinate methods, was not possible due to the limitations on the 
computer time available. An intuitive prerequisite on the grid 
size is that it should be less than the mean free path. We also 
note that the integration of the differential equations of the 
moment and discrete ordinate methods is a first order scheme. Hence 
the accuracy of the solution increases as the grid size is decreased. 
A possible gauge to test the 'reasonableness' of the grid size is 
to check and see if the mass flux along the slot is conserved. The 
percentage deviations in the mass flux calculated can be taken as 
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an index on the necessity for more grid points. A trapezoidal 
integration scheme is used to integrate pu across the slot 
section to give the total mass flux through the slot. 

A check of the moment solutions indicates that with NN 
between 3 and 5, the maximum deviation varies, respectively, from 
3 to 1%. This is true both for the free molecule and transition 
flow calculations. For free molecule flows by the discrete ordinate 
method, with either trapezoidal or cotangent integration, the maxi- 
mum deviation is 4 to 2% for NN between 3 and 6, respectively. 
However, in the transition regime for NN between 3 and 5, the 
maximum deviations were 9% to 4%. In fact, NN was kept = 4 for 
the transition regime calculations by the discrete ordinate method. 

Figure 13 shows a plot of Qf m /Qt ran vs. [Kn Q /(l+a)] 
[Qcont/Qfm^ 5 where Q is the mass flux through the slot scaled 
by p o a(2RT 0 )^, and the subscripts fm, cont and tran denote 
the free molecule, continuum and transition flow regimes, respec- 
tively. This Correlation is in the spirit of the work of Sherman 
(1963). However, a single curve for all values of the Kn Q could 
not be obtained. This was presumably because the mass flux through 
a finite length slot in the continuum limit was not known correctly. 
The Q cont used was obtained from the Poiseuille formula valid only 
for infinite length slots. 

Figure 14-1 shows the velocity profile (pu) across the 
slot section at the slot entrance, exit and midsection for i =12, 
Kn Q =5,1 and 0.5, and a = 0.1 . 

Comparing the moment and discrete ordinate methods, the 
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discrete ordinate method seems to give good results for all the 
values of parameters considered. However, the amount of computing 
involved is many times that by the moment method (see below). The 
free molecule solution by the moment method is deficient, but the 
transition flow solution by the moment solution quickly merges 
with that by the discrete ordinate method. For the same number of 
grid points NN and slot length, the computation time for each 
iteration by the moment method varies as £(NN) 2 , while by the 
discrete ordinate method as A (NN) . Further, it is found that 
in general the discrete ordinate method involves a greater number 
of initial value iterations since it is very difficult to guess the 
values of the g ( 1 ,1 ;£ a ,3ir/4) with any precision. 

As Kn Q is decreased, the number of grid points has to 
be increased since physically it seems that the distance between 
the grid points has to be less than the mean free path. Thus both 
our methods of computation are going to be more and more time con- 
suming (on the computer) as Kr\ n is decreased. Further, as Kn 
is decreased, the iteration for the initial values becomes time 
consuming, but we also found that the iterations tend to diverge 
very quickly if the initial values given are not judiciously chosen. 
The only way to avoid divergence seems to be to decrease Kn Q 
slowly. Even this is not sufficient at times. The iterations have 
to be nursed along very carefully by decreasing e (see Appendix 
C) appropriately (e is the increment in the initial values). 

5.3 Improvements and Suggestions for Future Work 

As discussed earlier, the solutions for the wall flux, etc. 
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seem to have a periodicity for even integer values Of x. There 
is no reason to expect this in an actual flow. The reason it 
occurs in our solution is that tt/ 4 is the only oblique (with 
respect to the slot axis) discrete ordinate value of 6 chosen. 

The obvious remedy is to take successively more numbers of oblique 
discrete values of 0. Unfortunately, the computer time involved 
per iteration goes up very rapidly. 

In view of the computer time consumed by discrete ordinate 
method calculations, it seems more appropriate to use a refined 
moment method, one preferably in the spirit of the work of Lees 
(1959). We do know the distribution function at any point in the 
slot for free molecule flow. In the solid angles subtended at the 
point by the slot entrance and exit, f, the normalized distribu- 
tion function has the values exp(-£ ) and a exp(-£ ). In the 
solid angle subtended by each wall, since the wall flux varies 
linearly with x, the normalized distribution function f is 
given by f = (a^ + b^ cot 0)exp(-£ ), where a^ and b^ are 
known functions of x arid y. These distribution functions could 
be generalized to include unknown functions of the physical space. 
The advantage of this method is that it could represent the free 
molecule distribution very well. Also, in view of the good results 
of the discrete ordinate method with the cotangent integration, in 
the transition regime, it seems reasonable to expect favorable 
results towards the continuum regime. More important, the computer 
time required for computation would only be of the same order as 
that by the simple moment method. 
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6.0 CONCLUSIONS 

1. The moment method is used to obtain the free molecule and transi- 
tion flow solutions. The range of parameters considered is: 

a = 0.1, 0.5 and 0.8, Z = 1 , 2, 4, 8 and 12, and Kn Q = 0.5, 
1.0, 5.0, °°. The free molecule solution seems to be significant 
only in the range & = 6 to Z = 12. However, the errors in 
the calculated value of the mass flux compared to the solution 
of Reynolds and Richey (1967) varies from 25% for Z = 6 to 
35% for Z = 12. The errors in fact get progressively worse as 
the slot length Z is increased beyond z = 12. The transition 
flow characteristics are vey well portrayed. For Kn Q =0.5 and 
1.0, Z = 1, 2, 4, 8 and 12, and a = 0.1, 0.5 and 0.8, the 
calculated value of the mass flux agrees to within 5% of the 
value of the mass flux obtained by using the plane Poiseuille 
formula with slip boundary conditions. 

2. The discrete velocity ordinate method is used; the velocity 
field is represented by two speeds and eight discrete points 
in the angular space. The integration over the speeds of the 
velocity field is done by a Gaussian quadrature. The integration 
over the angular space is done by two methods; first by a simple 
trapezoidal integration, and then by a cotangent integration 
that accounts for the special discontinuous nature of the dis- 
tribution function in the free molecule limit. The free molecule 
and transition flow solutions are obtained by both forms of 
angular integration. The range of parameters considered is the 
same as for the moment method. Compared to the solutions of 
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Reynolds and Richey (1967), the free molecule solution by the 
discrete velocity ordinate method with trapezoidar integration 
overestimates the mass flux by as much as 25% for l = 12. 
However, the cotangent integration in view of its special con- 
struction represents the free molecule Solution very well. 

For the transition flow calculations, the results by the two 
methods of angular integration vary less and less as Kn Q 
decreases and l increases. For the values of a arid SL 


considered, the agreement with the moment method transition 

flow solutions is very good, having a maximum difference of 

6 % for Kn = 5.0. This difference decreases as Kn„ is 
o o 

decreased. The disadvantage of the method is the large amount 
' of computer time required for the calculations. 


3. For the range of parameters considered, the calculated value 
of the mass flux through the slot, m, can' be 'represented to 
within 5% by the following equation: 


m_ 

P 0 a(2R T 0 ) 


where. 


1/2 a 


_ 1 , 0.1331 


Kn. 


t i 


Kn- 


I XI i 

(l+a)T-| + 2(T 1 - -r 1 °g s i ) 


Kn * . / 0 2 , /j \ 1/2 

+ ^72 log S 1 log ( ) } , 


■r _ (1-a)^ 

1 1 " l + 4.55 - 2.85 exp{0.21(Kn Q -5)} 


2 Kn f 8 ( 1 +a+T ) 


1 2Kn Q + 8 (l+a-T) 
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■ -i 


R is the gas constant, 'a' the slot half width,; p Q and T Q 
the upstream density and temperature, respectively. 
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FIG (Q SLOT GEOMETRY 
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FIG- 2 GRID POINTS FOR THE MOMENT 
METHOD 
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FIG-3 GRID POINTS FOR THE DISCRETE 
ORDINATE METHOD 
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FIG-4 COMPARISON OF SLOT MASSFLUXvsl 

BY THE VARIOUS SOLUTIONS FOR Kn=00 
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FIG -5 FREE MOLECULE WALL FLUX vsX,l=4,cc=0 
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APPENDIX A 

EXISTENCE AND UNIQUENESS OF MOMENT SOLUTIONS 

The source material used here is Saranson's work (1962). 

Let (n) be a vector function of (xj ,x 2 x^). In a region 

R, let it satisfy the boundary differential equation. 


9n . 


,2, A ij ^ * c j = 0 


On the boundary S of R, let the boundary condition (H)(n) 
be specified. 

Let the following conditions be satisfied: 

1) The matrices (A..) are symmetric. 

J 

2) The matrices (A-^.) are nonsingular. 

3) and (C^) are bounded matrices. 

4) The roots of the characteristic equations are real. 

5) Define, 


(g) 


1 N ^ i 

00 = 1 ((C) + (C)' - J ) 

c i=l 3x i 


where the prime sign 1 is used to denote the transpose of 
the matrix. Then, j ( (k) + (k) ' ) > 0. (If this condition , 
is satisfied, the system is said to be a positive symmetric , 
system. ) 

6) Call the unit exterior normal to the surface S as 
(s) = ( 0 q »© i »• . . >0 r ). Along S, define the matrix 


< b > =4 i 

c i=0 1 1J 


9 
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and write, (B) = (B + ) + (B_) and (m) = (B + ) - (B_). Then 
(in) has a nonnegative symmetric part. 

7) If the above six conditions are satisfied and (N) = (m) - (B) 
with (g) in the range of (N), then (N)(n) = (g) is said 
to be an admissible boundary condition. That is, for the 
given differential equation and boundary condition, a unique 
solution exists. 


Considering our problem, define 


and 


/if 

Kn 0 (T/T 0 ) 


'00 


CT 


(M) = 


M 


'01 

10 


L M 11 



Neglect compared to in the last equation of the 
set of equations (3.11). Our equations and boundary conditions 
then have the form, 


+ 




0 

0 

1 

0 

0 

0 

0 

0 


0 

1 

0 

0 

0 

0 

0 

0 




(M) = 


0 

1 

0 

0 


0 

0 

1 

0 



9 


0 


(A. 1 ) 
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wi th , 



as the unknowns, where 



and m, n, p, q and r are arbitrary functions whose choice is 
at our discretion. 

Substituting for (M) in Eq. (A.1 ) , we have 
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C ( A-j ) 3 I+ (A 2 ) ^ + (0^3(6) = 0 , (A. 3) 

where 




= 8x + 3y + Q 




Condition 1) is satisfied since (A-|) and (A 2 ) are symmetric 

matrices. By an appropriate choice of Z, m, n, p, q and r we 

can make (A^ ) nonsingular and (A-j ) and (C-j ) bounded. 

Conditions 2) and 3) are thus satisfied. 

Calculation of the roots of the characteristic equation 

yields values of 1 and -1 --both of which are real. Thus 

condition 4) is satisfied. Further, we have ^ ((k)+(k)') = 

((C-j )+(C-| ) 1 ). Thus if we choose Z, m, n, p, q and r such 

that Z , Z , m , m r , r are all positive, then 

x y x y x y 

j ((k)+(k)') = ( (C-| )+( C-j ) 1 ) is greater than zero. Thus 

condition 5) is satisfied, i.e., the system is symmetric 
positive. 

Finally, choose Z, m, n, p, q and r such that they are all 
positive in our domain, and also such that p > 2/rr m > 4irr and 
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n > 2 /jt A > 4irq. Then conditions 6) and 7) are immediately 
satisfied. 

A choice of A, m, n, p, q and r can easily be established 
such that the above conditions are satisfied. Thus a unique 
solution exists for (G). But since (M) is a simple algebraic 
combination of (G), hence a unique solution exists for (M) as 
wel 1 . 
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APPENDIX B 

FREE MOLECULE FLOW RESULTS FOR THE MOMENT METHOD 


Let us first consider the case where i is an even integer, 
equal to 2n, where n is an integer. Let us define 


4(1+Jt/ir) 


Then, 


for 


while 


for 


H(S) 

= D//F 


for -1 < S < 2n+l 

G(S) 

1 

2 

- D 


for -1 < S < 1 

G(S) 

a 

2 

- D 


for Jl-1 < S < 2n+l 

G(S) 

1 

= 2 

- {1 

+ i 

7T 

(m+1 )} D 

l+2m 

< S < 

l+2(m+l) 

with m = 0,1 ,. . . ,( 

G(S) 

1 

2 

- {1 

+ i 

7T 

(n-l)> D 

2n-3 

< S < 

2n-l 




Secondly, let us consider the case where i is an odd integer, 
equal to (2n+l), where n is an integer. 

Let us first define the following quantities: 

(1-g) 


E = 


4(1 + — ) 

7T 1 
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and 


2ir(2n+3) + 4(2n 2 +4n+l) 

Then, 


H(S) 

= E/Zif 

for 

S = -1 ,0,1,... ,(2n+2) 

9 

H(S) 

= (n+1 ) /if B 

for 

2m > S > 2m-l 




with 

m = 0,1 ,2,... ,n+l 

9 

H(S) 

= (n+2) /if B 

for 

2m+l > S > 2m 




wi th 

m = 0,1 ,2,. . . ,n-l 

9 

while. 





G(S) 

II 

M— ' 

m 

for 

S = -1 ,0,1 

9 

G(S) 

- i - E - — E 

2 7T 

for 

S = 2m or 2m+l 




with 

m = 1 ,2,.. . ,n 

9 

G(S) 

- ^ J_ c 

" 7 + E 

for 

S = 2n+2 

9 


G(S) = j B - 4m(n+l ) B 


for 2m > S > 2m-l 
with m = 0,1,2,... ,n+l , 

and 

G ( s ) = \ \ {(2n+3)u+2}B - 4m(n+2)B 

for 2m+l > S > 2m 
with m = 0,1 . ,(n-l ) 

As noted in the text, solutions for G and H can be obtained 
for any rational value of i. Except for integer values of 
the expressions for G and H are very complicated. 
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APPENDIX C 

ITERATION SCHEME FOR INITIAL VALUES 

Using either the moment or discrete ordinate method we have a 
problem where, in the discrete formulation, we have to determine a 
simultaneous set of values q-| ,q 2 %i suc ^ that 

(p 1 ,p 2 ,...,p m ) = M[(q 1 ,q 2 q m )] 

or 

£ = M(g_) (C.l) 

The nonlinear operator M represents the marching procedure 
described in the text, £ represents the given boundary values 
at x = SL (and the wall for the discrete ordinate method) and 
£ the unknown boundary values at x = 0 (and the wall for the 
discrete ordinate method). 

To solve (C.l) we use an iterative method based on linear- 
izing Eq. (C.l). Let £ n be an approximation to the solution 
with 

£ n = M( a n ) (C.2) 

If £ n is close to the true solution £ it seems reasonable 
to assume that £ - £° is, to good approximation, linearly dependent 
on (£ - £ n ). We therefore calculate £ n+ ^ using 

FT1 

(P. n+1 - P. n ) = l A..n (q.n+1 - n) (C.3) 

i i j= l ’3 3 

The elements of the matrix A. . n are given for i = l,...,m 

* J 

by 
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n 



5 


(C.4) 


where Ap^ 11 is the change produced in p.. by varying q^ 11 to 
q. n + Aq. while holding all other members of £ n fixed. To 
determine all the members of A. . n we now vary j from 1 to m. 

* J 

We have to "march" through the slot m times to determine the 
matrix A. . n . 

J 

The solution of Eq. (C.3) to obtain q. , for j = 1 to m, 

J 

was done using the University of California Computer Center Library 
Subroutine LINEQF. An optimum value of Aq . (j = 1 to m) exists 

J 

for smooth and quick iteration to the initial values. All the 

Aq j 1 s (j = 1 to m) were set equal. The final choice of the single 

value, e, was found to depend very critically upon a and Kn Q , 

-4 -7 

ranging from 10 to 10 . A smaller value was required if a 

was increased and/or Kn Q was decreased. 
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APPENDIX D 

r 2 

GAUSSIAN QUADRATURE FORMULAS FOR / dc ce" c h^c) 

o 

A general reference for the material here is Kopal (1961). 

We consider quadrature formulas of the form 

/ w^c) ^(c) dc ~ l h 1 (c |< ^ n ^) , . (D.l) 

cl K— 1 

where and C|^ n ^ are the Gaussian quadrature weights and 

roots, respectively , in the interval a < c < b. 

Define an inner product by 
b 

(h|,g-|) / w^c) hi (c) g,(c) dc . (D.2) 

a 

There exists a set of polynomials 

P n (c) = A n c n + B n c n_1 + ... (A n f 0) (D.3) 

which are mutually orthogonal with respect to this inner product, 
that is , 

(Pi»Pj) =0 , for i f j . (D.4) 

These orthogonality conditions define the polynomials up to a 
multiplicative constant which for our purpose is set equal to 
unity. 

The Gaussian quadrature roots c are the roots of 
P n (c). The Gaussian quadrature weights are given by 
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(n) _ 


(p n-l * p n-l ^ 


TTTTnTy 

i [C K } 


W 


D (C 

p n-r c K ' 


(D.5) 


We are interested in evaluating an integral of form 
00 ~2 

/ dc (c e' c ) h 1 (c) , 

o 

where h.(c) is a polynomial in c of degree 3 or less. Trans- 

-c 2 

lated in the symbols involved in Eq. (D.l), w^(c) = c e and 
n = 2. Thus using Eq. (D.4), we get by successive calculations, 


% ■ 1 

pi - <=- 4 

p 2 ■ {c • 4 ( )} p i - (i ■ T > 


Thus, the required Gaussian quadrature roots are the roots of the 
equation. 


p 2 = 0 

Solving, we get 


( 2 ) _ ( 2 ) _ 


1/2 

TT 

2T4^ffy 


(67r 2 -397r+64) 1/2 

2[4^j 


* 1.5181767, 0.54664006 


Further, using (D.5), we get the following values for the 
Gaussian quadrature weights, 
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i (4 -tt) 3 


u (2) . „ 

1 (6tt 2 -397t+64) 1 ^ 2 ( (6tt 2 -39tt+64 ) 1 /2 -tt ] /2 (tt-3) ) 


= 0.32523208 


and 


( 2 ) 


l (4-tt)' 


(6Tr 2 -397r+64)^ 2 ( (6 tt 2 -39tt+64 ) ^ 2 +tt^ ^ 2 (tt-3 ) ) 


= 0.17476792 


In terms of the notation used in the text, w a = H^ 2 ^, 


«b " H 2 (2) - 5 a = c/ 2) and 5 b = c 2 (2) . 



